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Abstract: Inverse source reconstruction is the most challenging aspect of 
bioluminescence tomography (BLT) because of its ill-posedness. Although 
many efforts have been devoted to this problem, so far, there is no generally 
accepted method. Due to the ill-posedness property of the BLT inverse 
problem, the regularization method plays an important role in the inverse 
reconstruction. In this paper, six reconstruction algorithms based on l p 
regularization are surveyed. The effects of the permissible source region, 
measurement noise, optical properties, tissue specificity and source 
locations on the performance of the reconstruction algorithms are 
investigated using a series of single source experiments. In order to further 
inspect the performance of the reconstruction algorithms, we present the 
double sources and the in vivo mouse experiments to study their resolution 
ability and potential for a practical heterogeneous mouse experiment. It is 
hoped to provide useful guidance on algorithm development and application 
in the related fields. 
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1. Introduction 

In the past decade, bioluminescence imaging (BLI) has been one of the hot topics in optical 
imaging and has been successfully applied in tumorigenesis studies, cancer diagnosis, 
metastasis detection, drug development, gene therapy, etc. [1—5]. However, the intrinsic 
planar imaging mode of BLI is just a qualitative imaging tool, which provides no absolute 
quantitative information and detects limited depths of a few millimeters inside tissues [4,5]. 
The bioluminescence tomography (BLT), the three-dimensional (3D) tomographic counterpart 
of BLI, has received much attention because it can provide localization and quantitative 
information of the bioluminescent source inside a living small animal [6-9]. 

In BLT, the tumor cells of the living mouse model are usually transfected with a luc-gene. 
Then, the oxidation of the substrate and luciferin results in a red-shifted light emission with a 
wavelength in the range from 500 nm to about 760 nm [3]. Based on the surface 
bioluminescent images and the heterogeneous tissue structure information captured by, for 
example, the micro-CT imaging system, BLT can reconstruct the location and distribution of 
the internal embedded bioluminescent light source and acquire quantitative volumetric 
reconstruction information [8]. 

The inverse source reconstruction is crucial for recovering three dimensional information 
of the light source. Many efforts have been devoted to studies on efficient inverse 
reconstruction [10-28]. Generally, some prior knowledge, such as the permissible source 
region and multispectral information, are used to reduce the ill-posedness of inverse 
reconstruction [10-14,16,17]. Meanwhile, the development of a fast and robust algorithm is 
another challenging problem. The common schemes include the regularization method, 
iterative method, stochastic method and so on [15,18-26]. Since the compressive sensing (CS) 
theory was introduced in 2006, which is a signal processing technique for efficiently acquiring 
and reconstructing images with sparse representation from far fewer measurements or signals 
below the Nyquist sampling theorem demands, it has been utilized by various areas of applied 
mathematics, computer science, electrical engineering and medicine [29,30]. As a new 
imaging technique, the sparse bioluminescent source and insufficient measuring signal are the 
characteristics of BLT, which can benefit by using the CS technique. Along with the 
development of the compressive sensing (CS) theory, the l v regularization method has become 
the mainstream strategy to obtain the optimal solution of the BLT inverse problem [17,22,24- 
26]. 

Generally speaking, regularization is used to construct the approximation of an ill-posed 
problem by a family of neighboring well-posed problems [31]. So far, the regularization 
methods have been used for solving various types of inverse problems in mathematics, 
statistics, machine learning, signal and image processing, heat conducting, and so on. The 
most commonly used form of regularization is the Tikhonov-type regularization. This method 
achieves a compromise between the minimization of the residual norm and the penalty term 
[32]. The / 2 -norm is the integrant penalty term for the Tikhonov-type regularization. The l 2 - 
norm penalty term imposes a small weight on small residuals, but a large weight on others. As 
a result, the / 2 -norm approximation has many modest residuals and relatively few larger ones, 
which generates smooth solutions and brings multipseudo sources to surround the true source 
in BLT reconstruction [22]. Considering the source sparsity characteristics and the 
insufficiency of the measured data, the /[-norm penalty term has been adopted in the 
regularization method and was successfully applied to BLT [24-26,33]. Unlike / 2 -norm 
regularization, / r norm regularization tends to produce a sparse solution with a large number 
of components equal to zero. However, in many practical applications, the solutions of the l\ 
regularizer are often less sparse than those of the / 0 regularizes The l 0 regularizer aims to find 
the number of non-zero elements in the solving domain, but it is an impossible mission in 
practical applications [34]. To find more sparse solutions than the l\ regularizer, / p (0 < p <1) 
regularization methods have been developed and employed to solve the minimization 
problems [34-36] and the BLT reconstruction problem [28]. 
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Although the / p (0 < p < 1, p = 2) regularization methods have been applied to the BLT 
reconstruction problem, so far, there is no generally accepted method which can be suitable 
for all of the reconstruction cases. This paper intends to fill the gap in the existing studies to 
systematically benchmark the performance of the / p -regularization-based BLT reconstruction 
algorithms. Our group has been conducting research in the BLT reconstruction method for 
years. In this paper, we evaluated six reconstruction algorithms, including the Truncated Total 
Least Square method (TTLS) [27], the Incomplete Variables Truncated Conjugate Gradient 
method (IVTCG) [25], the Truncated Newton Interior-Point method (TNIPM) [24], the 
Primal-Dual Interior-Point method (PDIP) [26], and the Weighted Iterative 
Shrinkage/Thresholding algorithm (WISTA) [28] that were developed by our group, and the 
Tikhonov regularization that we used for the BLT reconstruction in its standard way [31]. In 
order to investigate the response of these algorithms to the permissible source region, 
measurement noise, optical properties, and tissue specificity, we conducted a series of single 
source numerical phantom experiments. Then, the double source numerical phantom 
experiment and the in vivo mouse experiment were carried out to further test their 
performance. 

This paper is organized as follows. The next section describes the diffusion approximation 
model of BLT and the inverse reconstruction formula. The essence of the six regularization 
algorithms are then introduced briefly. In Section 3, the response of these algorithms to 
several influence factors are investigated and discussed based on the numerical phantom and 
in vivo mouse experiments to demonstrate the performance of the different algorithms. 
Finally, the discussion and conclusions are addressed in Section 4. 

2. Methodology 

2.1. Inverse reconstruction formula of BLT 

By assuming light as photons, the propagation of light in biological tissues is an extremely 
complex process that involves a series of optical behaviors, including the absorption, 
scattering, internal reflection and transmission. The radiative transfer equation (RTE) is often 
used to analytically model the process of light propagation. Although the accuracy of the RTE 
model is indubitable, it is generally difficult to directly solve it in complex biological tissues. 
Furthermore, the complicated implementation in the numerical settings and the large numbers 
of the resulting equations makes its efficiency low for BLT [37]. In infrared optical imaging, 
most of the biological tissues exhibit the characteristics of high scattering and low absorption. 
As the first order approximation model of RTE, the diffusion approximation (DA) is widely 
used to depict the process [8]: 

-V-(5(r)V(D(r)) + A (r)<D(r) = 5(r) (reQ) (1) 
<t(r) + 24 1 (r)^(r)(v(r)-V<l)(r)) = 0 (r e 5Q) (2) 

where S(r) is the source distribution, r is the position vector, <t(r) is the photon flux 

density, // a (r) is the absorption coefficient, Z)(r) = l/(3(// a (r) + (l-g)// s (r))) is the diffusion 

coefficient, jj s (r) is the scattering coefficient, g is the anisotropy parameter, and Q. cz R 1 is 

the domain of interest. Equation (2) is the Robin boundary condition, and v(r) is the unit 

outer normal to dQ. at r, A n (r) - (1 + R(r))/(l- R(v)) is the boundary mismatch factor and 

R(r) depends on the refractive index n of the surrounding medium that can be approximated 

by £«-1.4399>T 2 +0.7099«~ 1 +0.6681 + 0.0636« [38]. In the optical imaging experiment, 
the outgoing flux density on is measured by the highly sensitive CCD camera and can be 
expressed as follows [8]: 
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g(r) = -fl(r)(v(r).VO(r)) = -_<D(r)(r s SQ) (3) 

To solve Eqs. (1) and (2), the finite element method (FEM) was adopted and the equations 
were discretized by the finite element basis functions. A linear relationship between the 
measurable boundary data b and the unknown source distribution S is then established and 
the BLT forward problem can be obtained as [39]: 

AS = b (4) 

where the weight matrix A establishes the relationship between the measurements and the 
unknowns. 

Because the measurement is insufficient compared with the whole solving domain, the 
inverse reconstruction is an ill-posed problem. As a result, it is difficult and impossible to 
solve Eq. (4) directly. Usually, the regularization method is employed to convert Eq. (4) into 
the following optimization problem and S is replaced by x to make the notation more 
general: 

min ||Ax - bf + X llxll (5) 

II 112 II lip v ' 

n 

where ||x|| =(^]|x,.| P ) 1/p ar, d P denotes the parameter used in the l p regularization and satisfies 

0 < p < 1 or p = 2. Selecting an appropriate optimization method, the localization and 
distribution of the bioluminescent probes can be obtained from Eq. (5). 

2.2. Regularization methods 

2.2.1. Tikhonov regularization 

Tikhonov regularization method is one of the most commonly used methods to solve the ill- 
posed inverse problems. When p = 2, Eq. (5) is converted into the following minimization 
problem: 

min I Ax - bf 2 + X \\xf 2 (6) 

where X>0 is a properly chosen regularization parameter and controls the weight given to the 
minimization of the regularization penalty term relative to the minimization of the residual 
norm term. In the BLT application, the right-hand side b is always contaminated by various 
types of errors, such as the measurement errors, approximation errors, and rounding errors. 
Assuming the right-hand side b satisfies the discrete Picard condition, an explicit solution of 
(6) is given by 

x = (A 7 'A + Jiiy A T b (7) 

By using the singular value decomposition (SVD) explicitly, the system matrix A can be 
expressed as A = and the Tikhonov regularized solution is given by 

"afufb 
i=i cr,. + A cr. 

where U = (u l ,....,u n ) and V = (v, vj are matrices with orthonormal columns, 

U T U = V T V = I n , and where Z = diag{a x ,...,a n ) has non-negative diagonal elements 
appearing in non-increasing order such that rr, > • • • > a n > 0 [40]. 
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2.2.2. Truncated Total Least Squares Method 

The Tikhonov regularization method only considers the measurement errors on the right-hand 
side b , but the system errors also exist in the computed coefficient matrix of the system 
equation. In fact, whether the system errors come from the coefficient matrix or the 
measurement noise, they are inevitable in a real BLT experiment, so He et al. proposed the 
truncated total least squares (TTLS) method combined with a practical parameter-choice 
scheme termed as improved generalized cross-validation (IGCV) to cope with the BLT 
reconstruction problem including both the measurement noise and the system errors [27]. 

The TTLS method takes its basis in an SVD of the augmented matrix 

{A b)=UW T =%u t atf , <7 e R mx( " +1) , F e R ( ~ +1)x( " +1) , where U T U = V T V = I n+l , and 
2 = diag(o x ,...,o n ^ have singular values <r, >->a n+] >0 . The matrix V is partitioned 



such that V = 



( V V ^ 

' 11 ' 12 

V V 

V 21 *22y 



, V n e R" x * , and V 22 e R lx (" +1 *) , where k is the truncation 



V V T 

parameter. Then, the TTLS solution is given by x TTLS = -V l2 V 22 = — zr^r ■ Based on the 



l|2 

\\v 

22 2 



Mx— fell 2 



modified GCV (MGCV) criterion for TTLS G = — ^—r- , where enp 4 is the sum of the 

(m-enp k ) 

filter factors of TTLS, the truncation parameter k is obtained by the IGCV scheme. 

2.2.3. Incomplete Variables Truncated Conjugate Gradient (IVTCG) Method 

In view of the insufficient measurement, the source sparse characteristic and the 
oversmoothing of / 2 -norm regularization in BLT, linear system Eq. (4) is covert to the 
following / r norm regularization problem based on the CS theory: 

1 2 

min— +X be I (9) 

2 II 112 II 111 v > 

where ||x| = 'Y Jj x i is the l\ norm of x. Introducing vectors u and v, and the substitution 

x = u — v , Eq. (9) can be reformulated as the following convex quadratic program with 
nonnegative constrained conditions: 

mine 7 Z + — Z 1 Bz = F(z) 

2 (!0) 

s.t. z > 0 

where u>0 , v>0, z = [u vf , c = Al 2N + [-q qf , l 2N = [l,l,---,l] r , q = A T x, and 

B = [A - A] T [A -A]. Then, the incomplete variables truncated conjugate gradient (IVTCG) 
algorithm was presented by He et al. to solve Eq. (10) [25]. First, we define the working set 
T k by the violation of the Karush-Kuhn-Tucker (KKT) optimal conditions at the kth iteration, 

r k ={i\is {\,-,2N),[(z k \ > 0,(yF(z k )\ = 0]or[(z k ) ! = 0,(VF(z*)) ; . < o]} (11) 

where z k is the current iteration point. Then, we define the other two working sets I k and J k 
as follows: 
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/, = {/, e I k \ I < min {|/ t | , N, }} and J k ={j,eJ k \l< min {| J k \ , N max - N s }} (12) 

where 8 > 0 , A 7 , > 0 and N m!lx > N s are three constants, J k = T k \ I k and 

I k =|ze{l,---,27V}|[(z*). >0,(z*)./(VF(z*)). > <?]} . The variables with the set will be 

updated 7 t and J 4 at the current iteration. 

Finally, we need to decide the search direction. For the search direction <i* of I k , it can be 
obtained by solving the sub -problem 

rmn (VF(J)\ < + ^<X.< -^-«) (13) 
si.z* >0 

For J 4 , the search direction rf* = — w, . The other entries of d k are set to 0. After 

determining the search direction, the next iteration point is derived by the backtracking 
method. 

In the IVTCG method, the parameters N s and are set as = \^/^ J ar, d 



N„„ =N + 



. The time complexity of the IVTCG algorithm is 0(Mm ) + O(MN) . 



2.2.4. Truncated Newton interior-point method (TNIPM) 

Using the truncated Newton interior-point method (TNIPM), the object function (9) is 
transformed to a convex quadratic problem with inequality constraints: 

1 2 n 

min — ||Ax-Z>|L +1Ym. 

2 11 1,2 tr ' (i4) 

s.t. \x t \<Uj, i = !,-■■, n. 
Then, a logarithmic barrier function for the constraints |x,|<m, is added to the 
minimization function (14), and Eq. (14) is converted to the following problem: 

1 2 " 1 " 

min-||Ax-6|| 2 +^K,. — J][log(u i +x i )+\og{u r x i )} = F t {x,u) (15) 

2 ,=i t i= i 

where parameter t e (0, oo) . 

In order to compute the search direction, a preconditioned conjugate gradient method is 
adopted to compute the Newton system 



V 2 F,(x,u) 



Ax 
Au 



-VF,(x,u) (16) 



Then, we can obtain a feasible point by using the initial value and the search direction. 
The implementation of TNIPM in BLT reconstruction is presented in [24]. 

2.2.5. Primal-dual interior-point method (PDIP) 

As an efficient method for solving the optimization problems, interior-point methods have 
been applied in many areas in the past twenty years. In [26], the primal-dual interior-point 
method was adopted to solve the BLT inverse problem. Equation (4) is transformed to the 
following / r norm minimization problem: 
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mm |pt|| 
s.t. Ax-b 
x>0 



and then converted to the standard primal and dual forms in linear programming: 

Primal(P) : min c x Dual(D) : max b T y 

s.t. A T y + s = c 



(17) 



s.t. Ax = b 
x>0 



(18) 



s >0 



By introducing a logarithmic barrier term, (P) is converted to the following logarithmic 
barrier problems: 



P„ : min c T x-d^lrvcj 

j=i 

s.t. Ax = b 
x>0 



(19) 



where 0 is the barrier parameter. Then, we can obtain the Karush-Kuhn-Tucker conditions 



for P, : 



Ax = b,x > 0 
A T y + s = c 
1 



(20) 



e 



-XSe-e = 0 



where e = (l,...,l) r , the matrix X and S are the nxn diagonal matrix whose diagonal 
entries are precisely the components of x and s respectively. 

Hence, a new feasible point and the Newton direction (Ax, Ay, As) at the current value 

(x* , y k , s k ) can be obtained by solving the following equation system: 

AAx = b- Ax k 

< A T Ay + As = c-A T y" -s" (21) 
S k Ax + X"As = 0e-X k S k e 

The general framework of PD1PM is presented in [26]. 
2.2.6. Weighted iterative shrinkage/thresholding algorithm 

In some cases, the / r norm regularization methods are often less sparse than / p -norm (0 < p < 
1) regularization methods. In the paper [28], we proposed a weighted iterative 
shrinkage/thresholding algorithm (WISTA) to solve the / p -norm minimization problem (5) (0 
< p < 1). Because we cannot solve the / p -norm minimization problem directly, we reformed it 
by the weighted / r norm: 

min cox subject to Ax = b (22) 

A>0 

In order to obtain an ideal solution of (23), we used the iterative form as follows: 



: soft(x' +(A' (b - Ax' ' )) / a, coX I (2a)) 



(23) 
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where the weight vector a is calculated from the last iterate x' , soft(z, q) is the soft-threshold 

function with a threshold q, and soft(x, g)=sign(x)max(|x|-g, 0) . For the implementation of 

WISTA, we set two stopping conditions. The first is that the quantity of descending in the 
objective function is less than the stopping threshold. The second is that the algorithm reached 
the maximum number of iterations. 

The final algorithm of WISTA is formulated in [28]. 

3. Experiments and results 

In this section, a series of numerical experiments and an in vivo experiment were conducted to 
benchmark the performance of the six investigated reconstruction algorithms. We used two 
indicators to evaluate the reconstruction results: the reconstruction time (Time) of these 
algorithms and the distance (LocErr) between the actual source position and the 
reconstruction position. We defined the reconstruction position as the point with the 
maximum density of the reconstructed source. The regularization parameter of the Tikhonov 
method was chosen by the classical L-curve method and the parameter-choice scheme of the 
TTLS method was the improved generalized cross validation (IGCV). The regularization 
parameters of TNIPM and IVTCG used in reconstruction were manually optimized and they 
were set as le-6 in this paper. All experiments were performed in MATLAB on a personal 
computer with 2.66 GHz Intel Core 2 Quad CPU and 4 GB RAM. 

3.1. Numerical Experiments 

3.1.1. Reconstruction in a single source 

In this part, a cylindrical heterogeneous phantom was employed to demonstrate the features of 
the investigated algorithms. The numerical phantom consisted of muscle, lungs, heart, bone 
and liver, with the dimensions of 30 mm in height and 20 mm in diameter, as shown in Fig. 1 
[39]. The optical parameters of the five tissues were specified as listed in Table 1. A solid 
sphere source with a 1 mm radius was centered at (3, 5, 0) inside the right lung. Meanwhile, 
3874 nodes and 17763 tetrahedral elements were used in the inverse reconstruction. In order 
to compare the performance of the six algorithms among various conditions, we designed five 
experiments for the single source reconstruction. 




Heart Bone 



Muscle 



Source 




Liver 



Fig. 1 . Heterogeneous phantom with a single light source composed of muscle, lungs, heart, 
bone, liver and the source in the right lung. 



Table 1. Optical parameters of the heterogeneous phantom [39] 



Material 



Muscle 



Lungs 



Heart 



Bone 



Liver 



/4 [mm '] 
fj, [mm '] 
g 



0.010 
4.000 
0.900 



0.350 
23.000 
0.940 



0.200 
16.000 
0.850 



0.002 
20.000 
0.900 



0.035 
6.000 
0.900 
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1. Reconstruction using different permissible source regions 

As a priori knowledge, the permissible source region (PSR) is commonly used in the inverse 
reconstruction of BLT. Different PSR strategies may lead to different reconstruction results. 
In Fig. 2, we can see the surface flux distribution of the heterogeneous phantom in the 
coronal, sagittal and translucency views respectively and six red lines indicate the PSR. We 
set four different PSRs to verify its influence on the inverse reconstruction of the six 
algorithms. The PSRs were defined as follows: PS1 = {(x, y, z)| the left lung}; PS2 = {(x, y, 
z)| 0<x, y<7, -5< z <5}; PS3 = {(x, y, z)| -5<x, y<8, -5< z <5}; and PS4 = {(x, y, z)| the 
whole phantom}. The volume ratios of the four PSR were 3.26%, 5.20%, 17.93% and 100% 
respectively. The reconstruction results are listed in Figs. 36. Since TTLS could only resolve 
the overdetermined equation, it did not provide results for PS4. In Fig. 6(b), because the 
reconstruction time of PDIP in PS4 was too long so as to make the time curves almost 
indistinguishable in the Cartesian coordinate system, we used the log plot to depict the 
reconstruction time of various methods. 

We can draw the conclusion that the smaller the PSR, the better the reconstruction results 
obtained. From Figs. 3-5, we found that the reconstruction results of / p (0 < p < 1) 
regularization methods (TNIPM, IVTCG, PDIP, WISTA) were better than those of l 2 
regularization methods (TTLS, Tikhonov). From Fig. 3, the reconstruction results of 1% 
regularization methods were scattered around the true source. On the contrary, from Fig. 4 and 




Y X Y 



(a) (b) (c) 

Fig. 2. The surface flux distribution of the heterogeneous phantom from different views, (a)-(b) 
The coronal and sagittal views, respectively, (c) The translucency view of (a). 




(d) (e) (f) (g) 



Fig. 3. Axial views of the BLT reconstruction results of h regularization methods using 
different PSRs at z = 0mm. (a)-(c) Results of TTLS; (d)-(g) Results of Tikhonov. 
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Fig. 5, the reconstruction results of Z p (0 < p < 1) regularization methods were more accurate 
and the reconstruction sources were smaller than those of l 2 , except for the results of PDIP. 

From Fig. 6, the smallest location error was obtained by PDIP for PS1 and the second- 
smallest location error was obtained by IVTCG for PS2. As the permissible source region 
became larger, the location error of PDIP increased. The worst result of IVTCG was obtained 
for PS1 and it may be due to the fact that IVTCG couldn't develop its capacity of sparse 
reconstruction in a small PSR. The reconstruction results of TNIPM and WISTA remained 
stable with different PSRs. The reconstruction times of WISTA and IVTCG were smaller than 
the remaining algorithms and PDIP was the most time consuming algorithm. On average, the 
reconstruction locations of IVTCG, TNIPM and WISTA outperformed the other algorithms. 



n it i 



(a) 



(b) 



(c) 



(d) 



□□□□ 



(e) 



(0 



(g) 



(h) 



□ □□□ 



(i) 



(j) 



00 



(1) 



Fig. 4. Axial views of the BLT reconstruction results of /] regularization methods using 
different PSRs at z = 0mm. (a)-(d) Results of TNIPM; (e)-(h) Results of IVTCG; (i)-(l) Results 
of PDIP. 



nnnn 



(a) 



(b) 



(c) 



(d) 



Fig. 5. Axial views of the BLT reconstruction results of WISTA using different PSRs at z - 
0mm. 
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(b) 

Fig. 6. Performance metrics for the six algorithms using different PSRs. (a) The distance errors 
of the BLT reconstruction results; (b) The reconstruction time of various methods. 

2. Reconstruction at different measurement noise levels 

In this experiment, random noise at different levels was added to the synthetic data to evaluate 
the stability and robustness of the six algorithms using the permissible source region of PS1. 
The reconstruction results under different noise levels (0%-50%) of the Gaussian noise added 
to the synthetic data are shown in Fig. 7 and Table 2 [41]. The noise was generated by a 
MATLAB function randn and each algorithm was tested with a different noise sample. 
Tikhonov, TNIPM, 1VTCG and WISTA were stable at all noise levels. The location error of 
TTLS increased a little bit at the noise level of 10%, and then it was kept stable until the noise 
level of 40%. The reconstruction results of PDIP behaved similarly. At the noise level of 50%, 
both the TTLS and PDIP methods became unstable, so we ran them several times and just 
showed the best results in Fig. 8(a) and Fig. 8(e) respectively. In Fig. 7(b), TNIPM has a 
larger computation cost and IVTCG methods are fast but produce larger localization errors. 
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Noise Level 




(b) 

Fig. 7. Performance metrics of various algorithms at different noise levels, (a) The distance 
errors of the BLT reconstruction results; (b) The reconstruction time of various methods. 



Table 2. Reconstruction results at different measurement noise levels 



Noise level 


Method 


TTLS 


Tikhonov 


TNIPM 


IVTCG 


PDIP 


WISTA 


10% 


Loc Err(mm) 


1.843 


2.084 


1.550 


2.779 


1.234 


2.084 




Time(s) 


6.314 


0.984 


11.289 


0.955 


7.128 


2.273 


20% 


Loc Err(mm) 


1.843 


2.084 


1.550 


2.779 


1.311 


2.084 




Time(s) 


4.423 


0.961 


13.090 


0.915 


8.858 


3.059 


30% 


Loc Err(mm) 


1.843 


0.406 


1.550 


2.779 


1.311 


2.084 




Time(s) 


4.933 


2.084 


14.329 


1.3080 


5.7387 


2.566 


40% 


Loc Err(mm) 


2.431 


2.084 


1.550 


2.779 


1.311 


2.084 




Time(s) 


6.281 


2.702 


14.237 


1.077 


6.066 


3.022 


50% 


Loc Err(mm) 


1.258 


2.084 


1.550 


2.779 


1.536 


2.084 




Time(s) 


4.020 


0.746 


11.265 


0.875 


6.017 


1.834 
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Fig. 8. Axial views of the BLT reconstruction results of six regularization methods at the 50% 
measurement noise level at z = 0mm. (a)-(f) are the results of TTLS, Tikhonov, TNIPM, 
IVTCG, PDIP and WISTA respectively. 

3. Reconstruction using different optical parameters 

In order to investigate the influence of optical property on inverse reconstruction, four 
experiments with different deviations of absorption and scattering coefficients were conducted 
by adding 10% perturbation to the absorption and scattering coefficients of different organs. 
The corresponding results are shown in Table 3 and Fig. 9. Therein, the symbols were defined 
as: PI = {10% deviation of /u a and ju s I; P2 = {-10% deviation of /u a and ju s V, p 3 = {10% 
deviation of fl a and -10% deviation of // s }; and P4 = {-10% deviation of jU a and 10% 
deviation of fi a }. From these results, we found that the results of P4 were better than those of 
the others and the best reconstruction location was obtained by IVTCG in P4. The reason 
should be that the optical parameters exhibit the characteristics of high scattering and low 
absorption in P4, which conforms to the application scope of the diffusion approximation 
model. From Fig. 9, IVTCG achieved the smallest location error and WISTA had the least 
amount of computation time. The optical property changes had little impact on the inverse 
reconstruction of WISTA and TNIPM. The l p (0 < p < 1) regularization methods (TNIPM, 
IVTCG, PDIP, WISTA) still obtained better reconstruction sources than the l 2 regularization 
methods (TTLS, Tikhonov). 

Table 3. Reconstruction results using different optical parameters 



Case 


Method 


TTLS 


Tikhonov 


TNIPM 


IVTCG 


PDIP 


WISTA 


PI 


Loc Err(mm) 


0.782 


1.830 


1.753 


1.956 


2.005 


1.955 




Time(s) 


4.741 


0.804 


10.003 


0.408 


11.223 


0.031 


P2 


Loc Err(mm) 


2.689 


1.416 


1.818 


0.470 


2.005 


1.955 




Time(s) 


3.793 


0.897 


8.158 


2.571 


12.821 


0.106 


P3 


Loc Err(mm) 


2.452 


1.543 


1.753 


1.956 


2.652 


1.955 




Time(s) 


3.760 


0.977 


9.532 


0.228 


6.262 


0.046 


P4 


Loc Err(mm) 


2.345 


1.543 


1.995 


0.470 


1.485 


1.955 




Time(s) 


4.774 


1.022 


8.198 


0.205 


5.171 


0.113 
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(b) 

Fig. 9. Performance metrics for the six algorithms using different optical parameters, (a) The 
distance errors of the BLT reconstruction results; (b) The reconstruction time of various 
methods. 



4. Reconstruction in tissue specificity 

In this part, the influence of tissue specificity on the inverse reconstruction was observed in 
four experiments. In the first case of Nl, the absorption, scattering and anisotropy coefficients 
of all organs were set the same as that of muscle to simulate a 3D homogeneous model. In the 
cases of N2-N4, we simulated three heterogeneous models (see Fig. 10). In the case of N2, the 
absorption, scattering and anisotropy coefficients of the lungs, heart and liver were set the 
same as those of muscle, thus the model was divided into two parts: muscle and bone. In the 
case of N3, the absorption, scattering and anisotropy coefficients of the lungs and heart were 
set the same, that is, the model was divided into four parts. The case N4 was the PS1 in the 
PSR experiment and there were five tissues in the model. The reconstruction results are 
presented in Table 4 and Fig. 1 1 . From these results, we found that the reconstruction results 
of the case N4 were better than the other three cases for / p (0 < p < 1) regularization methods 




(a) (b) (c) (d) 



Fig. 10. Tissue specificity models (various colors are for various segmented tissues). 
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(TNIPM, PDIP, WISTA), but this could not be observed for the h regularization algorithms 
(TTLS, Tikhonov) and IVTCG. It may be that there were more optical parameters in tissues in 
N4, so it was more approachable in reality than the other three cases. As shown in Fig. 1 1(a), 
the PDIP's performance was the best in the four cases and the reconstruction results of 
WISTA ranked second. In Fig. 1 1(b), WISTA still used the minimal amount of time, followed 
by IVTCG and Tikhonov. 




fti N2 N3 N4 

(b) 

Fig. 1 1 . Performance metrics for various algorithms of different tissue specificity, (a) The 
location errors of BLT reconstruction; (b) The reconstruction time of various methods. 



Table 4. Reconstruction results for tissue specificity 



Case 


Method 


TTLS 


Tikhonov 


TNIPM 


IVTCG 


PDIP 


WISTA 


Nl 


Loc Err(mm) 


3.061 


2.084 


3.443 


2.680 


1.311 


1.955 




Time(s) 


3.687 


0.945 


12.664 


0.978 


5.063 


0.154 


N2 


Loc Err(mm) 


3.061 


2.084 


3.443 


2.680 


1.311 


1.955 




Time(s) 


4.802 


1.044 


12.859 


1.314 


5.167 


0.230 


N3 


Loc Err(mm) 


3.504 


2.084 


3.038 


2.084 


0.410 


1.955 




Time(s) 


5.413 


0.980 


21.104 


0.291 


6.348 


0.233 


N4 


Loc Err(mm) 


1.602 


2.084 


1.550 


2.779 


0.410 


2.084 




Time(s) 


14.881 


6.021 


14.592 


0.820 


6.883 


2.340 



5. Reconstruction at different source locations 

In this experiment, we carried out three experiments to observe the influence of the source 
locations on the inverse reconstruction. In case 1, the source location was (4,6,0) and the 
distance between the source center and the surface of phantom was 2.789mm. In case 2, that is 
PS1, the source location was (3,5,0) and the distance between the source center and the 
surface of phantom was 4.169mm. In case 3, the source location was (2,4,0) and the distance 
between the source center and the surface of phantom was 5.527mm. The reconstruction 
results are listed in Table 5 and the error bar can be seen in Fig. 12. From the results, we 
observed that the smaller the depth was, the better the reconstruction results obtained. The 
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standard deviation of WISTA was smallest indicating that it's robust to source locations. 
Except for PDIP, the l p (0 < p < 1) regularization methods (TNIPM, IVTCG, WISTA) 
obtained smaller standard deviations than l 2 regularization methods (TTLS, Tikhonov). 

Table 5. Reconstruction results at different source positions 



Case 


Method 


TTLS 


Tikhonov 


TNIPM 


IVTCG 


PDIP 


WISTA 


1 


Loc Err(mm) 


2.027 


1.593 


1.836 


1.582 


2.802 


1.582 




Time(s) 


4.004 


0.944 


10.686 


0.486 


10.687 


1.128 


2 


Loc Err(mm) 


1.602 


2.084 


1.550 


2.779 


0.410 


2.084 




Time(s) 


14.881 


6.021 


14.592 


0.820 


6.883 


2.340 


3 


Loc Err(mm) 


3.293 


1.089 


1.089 


2.941 


1.448 


2.272 




Time(s) 


5.443 


1.926 


10.980 


0.326 


7.362 


2.020 



Average Loc Krr(mm) 



O 
c 

W 2- 



O 



TTI.S Tikhonov TNIPM IVTCG PDIP WISTA 

Method 

Fig. 12. Error bar chart of the Loc Err at different source positions. 

3.1.2. Reconstruction of double sources 

In this part, three experiments were carried out to verify the performance of these algorithms 
on the double source resolution. In experiment 1, two sources were located in the left lung 
with their centers at (3, 5, 1.25) and (3, 5, -1.25), and with a distance of 2.5mm. In 
experiment 2, two sources were located in the left lung with their centers at (3, 5, 1.5) and (3, 




Fig. 13. BLT reconstruction results for l 2 regularization methods in a double source case in a 
3D view, (a)-(c) Results of TTLS; (d)-(f) Results of Tikhonov. 
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5, -1.5), and with a distance of 3mm. In experiment 3, two sources were located in the left 
lung with their centers at (3, 5, 2) and (3, 5, -2), and with a distance of 4mm. The right lung 
was specified as the PSR in the three experiments. The reconstruction results are presented in 
Figs. 13-15 and Table 6. From these results, we found that TTLS only reconstructs one source 
in case 1 and case 3, as shown in Fig. 13(a) and Fig. 13(c). In Fig. 13(f), the Tikhonov's 
reconstruction of the two sources was almost overlapping. On the contrary, l p (0 < p < 1) 
regularization methods (TNIPM, IVTCG, PDIP, WISTA) obtained much better results. 
Comparing Fig. 14 with Fig. 15, we found that the reconstruction sources of WISTA were 
smaller than that of the three /, regularization methods (TNIPM, IVTCG, PDIP), which 
confirm that WISTA was sparser than the /] regularization methods. 




Fig. 14. BLT reconstruction results of l\ regularization methods in a double source case in a 3D 
view, (a)-(c) Results of TNIPM; (d)-(f) Results of IVTCG; (g)-(i) Results of PDIP. 




Fig. 15. BLT reconstruction results of / p (0 < p < 1 ) regularization methods in a double source 
case in a 3D view, (a)-(c) Results of WISTA. 
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Table 6. Reconstruction results in double sources 



Case 


Method 


Source 


TTLS 


Tikhonov 


TNIPM 


IVTCG 


PDIP 


WISTA 


1 


Loc Err(mm) 


Source 1 


2.125 


0.603 


0.603 


1.631 


1.687 


1.851 






Source2 




1.719 


0.514 


1.166 


1.376 


1.746 




Time(s) 




3.623 


1.127 


6.237 


0.248 


3.075 


0.082 


2 


Loc Eir(mm) 


Source 1 


1.996 


1.988 


0.772 


0.772 


1.230 


1.777 






Source2 


2.029 


2.173 


0.626 


1.401 


1.493 


1.886 




Time(s) 




4.707 


1.189 


5.542 


0.262 


2.708 


0.125 


3 


Loc Err(mm) 


Source 1 






1.199 


1.199 


0.6262 


1.953 






Source2 


1.910 


1.369 


1.230 


1.010 


1.811 


1.634 




Time(s) 




6.142 


0.779 


13.822 


0.747 


2.043 


0.352 



3.2. In vivo mouse experiment 

The in vivo experiment was presented here to further verify the performance of these 
algorithms for the application of living animal studies. In the mouse experiment, an implanted 
source filled with luminescent liquid was placed into the abdomen of the living mouse. The 
source was 2.7 mm in diameter and 5.1 mm in length. The experimental details and the optical 
parameters of each organ were presented in [26]. In this experiment, the whole mouse was 
specified as the PSR, which made the weight matrix A in Eq. (5) undetermined. The 
reconstructed results are presented in Fig. 16 and Table 7. The TTLS did not provide results 
because it could not resolve the undetermined equation. From Fig. 16(c), we found that the 
reconstruction location of TNIPM was on the bottom left corner of the mouse model and we 
concluded that its reconstruction failed. From Table 7, we can see that WISTA obtained the 
best reconstruction results and the run time of PDIP was the longest. 

Table 7. Reconstruction results of the in vivo mouse experiment 



Method 


Tikhonov 


TNIPM 


IVTCG 


PDIP 


WISTA 


Loc Err(mm) 


5.004 


19.023 


5.004 


5.624 


4.535 


Time(s) 


57.597 


91.537 


444.986 


1241.296 


101.235 



* 




(a) 




Fig. 16. BLT reconstruction results of the in vivo mouse experiment, (a) the 3-D view of the 
segmented micro-CT slices of the imaged mouse with a luminescent source; (b)-(f) the 
reconstruction results of Tikhonov, TNIPM, IVTCG, PDIP and WISTA. 



4. Discussion and Conclusions 

In this paper, we provided a comprehensive assessment of six regularization algorithms 
developed by our group. In order to investigate the performance of these algorithms in the 
BLT inverse reconstruction, we carried out five single source experiments, a double source 
experiment and an in vivo experiment. Because the inverse reconstruction equation is not well 
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posed, the reconstruction results were influenced by various factors, such as the permissible 
source region, measurement noise, optical properties, tissue specificity and source positions. 
In this paper, we observed all of the five influential factors which affected the results of BLT 
reconstruction in the single source experiments. In the double source case, the ability of the 
six algorithms in resolving double sources was investigated. Finally, the feasibility and 
practicability of these algorithms were further validated by the in vivo experiment. 

The extensive experiments showed that, under a wide range of conditions, there isn't 
always a winner that achieves the best performance for all of the investigated cases. Various 
results for the six reconstruction algorithms under different conditions were obtained and 
some interesting conclusions could be addressed. Consistent with our anticipation, the 
reconstruction results became better as the PSR decreased, especially for the l 2 regularization 
methods (TTLS, Tikhonov), which could be mainly attributed to the fact that the weight 
matrix was overdetermined and the sparsity characteristics of the inverse problem didn't pan 
out. In the noise experiment, four algorithms (Tikhonov, TNIPM, IVTCG, WISTA) were 
robust to noise. The PDIP became unstable at a noise level of 50% and TTLS became unstable 
at a noise level of 40%. The optical property variation experiments suggested that we could 
get better reconstruction results if the mathematical model used to describe the light 
propagation in the biological tissue fits the optical characteristics of the tissues. The tissue 
specificity experiments demonstrated the advantage of the heterogeneous model over the 
homogeneous model, which has been verified by many studies. But for l 2 regularization 
methods (TTLS, Tikhonov), especially the Tikhonov method, the reconstruction results didn't 
improve significantly, so we should use the / p (0 < p < 1) regularization algorithms in the 
reconstruction of the heterogeneous model. In the last single source experiments, WIS AT was 
the most robust algorithm at different source positions than the other five algorithms. 

In the reconstruction of double sources, TTLS had two cases and Tikhonov had one case 
where they couldn't separate the two sources, which contrasted sharply with the 
reconstruction results of the l p (0 < p < 1) regularization algorithms (TNIPM, IVTCG, PDIP, 
WISTA). It is due to the fact that l 2 regularization is easy in generating smooth solutions, 
which leads to multi-pseudo sources that can't distinguish the two sources. 

In the in vivo mouse experiment, except for the TNIPM method, other algorithms 
performed normally. As in the double source experiment, WISTA still obtained the best 
performance because it could maintain its preferable sparsity characteristics. We saw that the 
reconstruction results of the in vivo mouse experiment were unsatisfactory compared to the 
simulation experiment. Actually, many factors influenced the reconstruction results of the in 
vivo experiment, such as the signal acquisition, image registration of the CT data and optical 
data, the quality of mesh dissection, the accuracy of the optical propagation model, the inverse 
reconstruction algorithm and so on. 

For most experiments, the / p (0 < p < 1) regularization algorithms (TNIPM, IVTCG, PDIP, 
WISTA) achieved better reconstruction results than l 2 regularization (TTLS and Tikhonov). In 
the four l p (0 < p < 1) regularization algorithms, IVTCG and WISTA were generally more 
effective than the other two methods. 

In order to evaluate the / p -regularization-based BLT reconstruction algorithms, we limited 
the review to six algorithms (TTLS, Tikhonov, TNIPM, IVTCG, PDIP, WISTA) that were 
developed by our group. This is mainly because our group has been engaged in the research of 
bioluminescence tomography for years, we have the source code of the six algorithms and we 
can show their best performance and give them a fair comparison. Some interesting 
conclusions were obtained in this paper which hope to provide useful information for the 
researcher in selecting and designing algorithms for BLT reconstruction. 
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